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EXECUTIVE  SUMMARY 


Numerical  simulations  of  the  deformations  and  breakup  of  drops  have  been  done.  Inertial  and 
viscous  effects  for  both  the  drop  and  the  ambient  gas  as  well  as  surface  tension  effects  are  fully  ac¬ 
counted  for.  The  simulations  help  determine  where  in  parameter  space  the  various  breakup  modes 
take  place,  how  long  breakup  takes,  and  what  the  resulting  drop  size  distribution  is.  The  goal  of 
the  investigation  is  to  provide  results  that  extend  and  complement  experimental  investigations,  and 
lead  to  better  engineering  models  of  drops  in  sprays.  Simulations  of  the  primary  breakup  of  a  flat 
interface  have  also  been  done. 

To  examine  the  breakup  of  drops  as  the  density  difference  becomes  small,  extensive  axisym- 
metric  simulation  of  four  systems  have  been  done:  Impulsive  and  gradual  disturbances  for  two 
different  density  ratios.  At  low  density  ratios,  the  density  disappears  as  an  independent  control 
parameter  and  it  can  be  shown  that  the  low  density  results  apply  to  density  ratios  as  high  as  two 
if  time  is  re-scaled  using  the  Boussinesq  approximation.  In  addition  to  full  simulations  where  the 
Navier-Stokes  equations  are  solved,  a  few  inviscid  simulations  have  also  been  done  to  isolate  the 
effect  of  viscosity. 

Four  nondimensional  numbers  govern  the  breakup  of  drops.  In  addition  to  the  density  and 
the  viscosity  ratio,  the  ratio  of  inertia  to  surface  tension  is  described  by  an  Eotvos  number  for 
gradual  disturbances  and  a  Weber  number  for  impulsive  acceleration.  The  effect  of  viscosity  is 
described  by  the  Ohnesorge  number  (the  ratio  of  the  viscous  force  to  the  surface  tension).  The 
simulations  have  resulted  in  a  fairly  complete  picture  of  the  evolution  at  small  density  ratios. 
For  small  Eotvos  and  Weber  numbers  the  drops  remain  spherical  in  all  cases,  independent  of 
the  Ohnesorge  number  and  the  density  and  the  viscosity  ratio.  If  the  Ohnesorge  number  is  low, 
the  deformations  of  the  drop  depend  only  on  the  Eotvos/Weber  number  (and  the  density  ratio). 
As  the  Eotvos/Weber  number  is  increased,  the  drops  deform  into  a  disk-like  shape  due  to  high 
pressure  at  the  fore  and  aft  stagnation  points  and  low  pressure  around  the  equator.  Increasing  the 
Eotvos/Weber  number  further  results  in  a  continuing  deformation  where  most  of  the  drop  fluid 
ends  up  in  a  torus  connected  by  a  thin  film.  For  moderate  Eotvos/Weber  numbers,  the  initial 
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momentum  of  the  drops  is  relatively  low  and  once  the  torus  is  formed,  the  rim  moves  faster  than 
the  film  for  drops  with  gradual  disturbances.  The  film  “bulges”  back  and  experimentally  it  is 
seen  that  this  bag  eventually  breaks.  The  simulations  have  shown  that  the  bag  break-up  mode 
is  a  viscous  phenomenon,  due  to  flow  separation  at  the  rim  of  the  drops  and  the  formation  of 
a  wake,  and  it  is  therefore  not  seen  in  inviscid  computation.  For  drops  subject  to  an  impulsive 
acceleration,  the  formation  of  a  backward  facing  bag  is  only  seen  for  the  higher  density  ratios.  Bag 
breakup  requires  a  driving  force  that  acts  stronger  on  the  drop  than  on  the  surrounding  fluid  and 
for  impulsively  accelerated  drops  this  driving  force  is  the  fluid  inertia.  As  the  density  difference 
becomes  small,  the  difference  between  the  drop  and  the  fluid  inertia  vanishes  and  the  low-density 
ratio  drops  simply  stop  and  surface  tension  pulls  them  back  into  a  spherical  shape.  Experimentally, 
bag  breakup  is  commonly  observed  for  impulsive  acceleration,  but  the  density  ratio  is  much  larger. 
Increasing  the  Eotvos  or  the  Weber  further  results  in  a  different  mode  of  breakup  that  also  depends 
on  the  density  ratio.  For  low  density  ratios,  the  fluid  initially  still  ends  up  in  the  rim  of  the  drop, 
but  the  initial  momentum  is  now  sufficiently  large  so  the  ambient  fluid  moves  the  film  faster  than 
the  torus,  leading  to  a  bag  that  extends  forward  of  the  drop.  For  higher  density  ratios,  not  all 
the  fluid  moves  to  the  rim  and  a  toms  that  is  connected  to  the  rest  of  the  drop  by  a  thin  sheet 
is  formed.  As  the  driving  force  is  increased  the  size  of  the  rim  is  reduced  and  for  very  high 
Eotvos/Weber  numbers,  small  drops  are  pulled  from  the  rim.  Examination  of  the  dynamic  of  the 
vorticity  generated  at  the  drop  suggests  that  the  shear  breakup  mode  where  fluid  is  stripped  from 
the  rim  of  the  drop  is  an  essentially  inviscid  effect.  In  the  transition  between  a  bag  breakup  mode 
and  shear  breakup,  drops  that  oscillate  in  a  chaotic  manner  are  sometimes  seen.  Such  transition 
phenomena  have  been  seen  experimentally  for  higher  density  ratios.  In  addition  to  the  Ohnesorge 
number  effect,  where  the  boundary  between  breakup  modes  is  shifted  to  higher  Eotvos/Weber 
numbers  as  viscous  effects  become  more  important,  the  fluid  and  drop  viscosity  can  change  the 
drop  shape  during  breakup  if  the  Ohnesorge  number  is  high  enough.  High  viscosities  can,  for 
example  lead  to  skirted  drops  at  low  density  ratios,  where  thin  fluid  skirts  are  pulled  from  the  rim 
of  the  fluid,  in  a  way  similar  to  the  shear  breakup  seen  for  higher  density  ratios. 

The  simulations  have  been  used  to  generate  ’’break-up”  maps  and  it  is  found  that  the  general 


11 


character  of  those  maps  agrees  with  what  has  been  found  experimentally  for  larger  density  ratios. 
At  low  Ohnesorge  numbers  the  transition  between  the  various  modes  depends  only  on  the  density 
ratio  and  the  Eotvos  number,  but  a  high  Ohnesorge  number  will  move  the  transition  to  a  higher 
Eotvos  number. 

Simulations  of  three-dimensional  aspects  of  the  drop  breakup  and  heat  transfer  have  recently 
been  initiated. 
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I.  BACKGROUND 


In  spray  combustion,  liquid  atomization  is  a  two  stage  process:  In  the  primary  breakup,  a 
liquid  jet  emerging  from  the  injector  breaks  up  into  drops  which  subsequently  undergo  secondary 
breakup  into  even  smaller  drops.  Secondary  breakup  increases  the  total  surface  area  of  the  fuel-air 
interface,  thus  enhancing  the  rate  at  which  the  fuel  evaporates  and  bums. 

Current  computational  models  used  for  engineering  predictions  of  spray  combustion  do  not 
resolve  the  motion  of  individual  drops.  Instead,  the  effect  of  the  drops  is  accounted  for  by  subgrid 
models,  computed  in  either  an  Eulerian  or  a  Lagrangian  way.  For  recent  descriptions  and  reviews, 
see  Drew  and  Passman1  and  Crowe,  Sommerfeld,  and  Tsuji2.  In  Lagrangian  models,  the  drops 
are  represented  by  point  particles,  that  can  be  split  into  two  or  more  particles  to  represent  drop 
breakup.  For  a  description  and  application  of  breakup  models  in  spray  combustion  simulations, 
see  Reitz  and  Diwakar3,  O’Rourke  and  Amsden4,  Liu,  Mather,  and  Reitz5,  Liu  and  Reitz6,  Kim 
and  Wang7,  and  Kong,  Han,  and  Reitz8.  Two  different  approaches  are  typically  used  to  model 
the  breakup.  The  Taylor  analogy  breakup  (TAB)  model  of  O’Rourke  and  Amsden4  is  based  on 
an  analogy  between  an  oscillating  and  distorting  liquid  drop  and  a  spring-mass  system  suggested 
by  Taylor9.  The  Reitz  wave  instability  model6,  on  the  other  hand,  is  based  on  a  linear  stability 
analysis  for  liquid  jets.  Both  of  these  simplified  models  contain  adjustable  parameters  that  must 
be  determined  by  experimental  studies. 

In  experiments,  the  drops  are  usually  accelerated  by  a  shock  wave  causing  a  step  change  in 
the  velocity  of  the  drop  relative  to  the  surrounding  fluid,  or  by  a  constant  body  force  such  as 
gravity.  The  results  are  generally  presented  in  terms  of  four  non-dimensional  parameters:  the 
relative  strength  of  inertia  and  surface  tension  which  is  characterized  by  the  Weber  number  for  an 
impulsive  acceleration  and  the  Eotvos  number  for  an  acceleration  by  a  constant  body  force;  the 
ratio  of  viscous  stresses  and  surface  tension  given  by  the  Ohnesorge  number;  the  density  ratio;  and 
the  viscosity  ratio  of  the  drop  and  surrounding  fluids. 

Early  experimental  studies  of  drop  breakup  due  to  impulsive  acceleration  include  those  of 
Lane10 ,  who  studied  the  shattering  of  liquid  drops  in  steady  or  transient  streams  of  air,  and  Hinze11 . 
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The  findings  of  Lane10  and  Hinze11  have  been  extended  to  a  broader  range  of  parameters  by 
Haas12,  Hanson,  Domich,  and  Adams13,  Ranger  and  Nicholls14,  Gel’fand,  Gubin,  Kogarko,  and 
Komar15,  Borisov,  Gel’fand,  Natanzon,  and  Kossov16,  and  others.  Krzeczkowski17  showed  that 
the  effect  of  drop  viscosity  is  not  significant  when  the  Ohnesorge  number  based  on  the  drop  prop¬ 
erties  is  less  than  about  0.1.  Pilch  and  Erdman18  examined  the  fragment  size  distribution  for  the 
so-called  bag  breakup  mode  and  found  that  it  was  made  up  of  a  large  number  of  small  fragments 
produced  from  the  burst  of  the  bag,  and  a  few  large  fragments  originating  from  the  annular  rim. 
Wierzba19  reviewed  the  literature  and  found  that  there  was  a  large  variation  in  the  reported  values 
of  the  critical  Weber  numbers  for  the  onset  of  the  bag  type  breakup.  His  own  experiments  showed 
that  small  changes  in  the  experimental  conditions  could  affect  the  drop  breakup  significantly.  Ex¬ 
perimental  investigations  of  the  breakup  of  falling  drops  have  typically  been  motivated  by  interest 
in  the  evolution  of  raindrops.  For  early  experiments,  see  the  study  of  the  breakup  of  large  drops  by 
Magarvey  and  Taylor20,  for  example.  The  secondary  breakup  of  liquid  drops  due  to  both  impul¬ 
sive  and  continuous  disturbances  has  been  examined  extensively  by  Hsiang  and  Faeth21  23  using  a 
shock  tube  and  drop  towers.  The  majority  of  the  data  are  for  atmospheric  conditions  (pd/ p0  >  500, 
Re  >  100,  Ohd  <  0.1),  although  a  limited  number  of  studies  for  smaller  density  ratios  and  higher 
viscosity  were  also  done.  For  the  impulsive  disturbance  case,  droplet  deformation  and  breakup 
maps  similar  to  those  produced  by  Hinze11  and  Krzeczkowski17  were  constructed  for  a  wide  range 
of  parameters.  Joseph,  Belanger,  and  Beavers24  studied  the  breakup  of  both  Newtonian  and  non- 
Newtonian  drops  in  a  high-speed  air  stream.  Their  experiments,  using  a  shock  tube,  resulted  in  a 
very  high  initial  acceleration  of  drops  and  the  authors  stated  that  the  Rayleigh-Taylor  instability 
was  the  primary  cause  of  breakup.  For  a  more  extensive  review  of  experimental  studies  of  sec¬ 
ondary  breakup  of  drops,  see  Clift,  Grace  and  Weber25,  Lefebvre26,  Bayvel  and  Orzechowski27, 
and  Sadhal,  Ayyaswamy,  and  Chung28. 

Most  of  the  experimental  studies  mentioned  previously  are  concerned  with  the  breakup  of  liq¬ 
uid  drops  in  air  due  to  impulsive  accelerations.  The  density  and  viscosity  ratios  are  much  higher 
than  those  considered  in  the  present  study.  While  those  experimental  results  are  not  directly  com¬ 
parable  to  our  simulations,  the  major  breakup  modes  remain  similar.  We  therefore  summarize  the 
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major  results  of  experimental  studies  of  impulsively  accelerated  drops  here.  When  the  Ohnesorge 
number  is  small,  the  effects  of  drop  viscosity  can  be  neglected.  At  low  W e,  the  drops  deform 
but  do  not  break  up.  As  the  acceleration  increases  past  a  critical  value,  the  drops  become  pro¬ 
gressively  flatter  and  eventually  break  up.  As  the  Weber  number  is  increased,  four  well  defined 
breakup  modes  are  observed  (see,  for  example,  Nigmatulin29): 

1.  Vibrational  breakup  mode  where  the  drop  disintegrates  into  two  or  four  equal-sized  smaller 
drops. 

2.  Bag  breakup  mode  where  the  original  drop  deforms  into  a  torus-shaped  rim  spanned  by  a 
thin  fluid  film  that  ruptures  into  tiny  droplets,  followed  by  disintegration  of  the  rim  into 
larger  droplets. 

3.  Shear  breakup  mode  where  small  drops  are  continuously  stripped  off  the  rim  of  the  original 
drop. 

4.  Explosive  breakup  mode  where  strong  surface  waves  disintegrate  the  drop  in  a  violent  man¬ 
ner. 

This  categorization  and  terminology  are  somewhat  arbitrary,  and  other  variations  have  been 
suggested  by  other  researchers.  For  example,  mode  3  has  also  been  called  “stripping-type 
breakup.”18,19  For  low  viscosity  drops  where  the  transition  process  shows  no  significant  depen¬ 
dencies  on  Ohd,  the  critical  Weber  number  is  approximately  10  for  the  first  transition,  20-60  for 
the  second,  and  1000  for  the  third.  These  numbers  should  be  considered  only  as  a  rough  guide 
because  there  are  large  variations  in  the  critical  Weber  numbers  in  the  available  experimental  data 
due  to  different  test  conditions.  At  higher  Ohd,  the  We  required  for  the  onset  of  deformation  and 
breakup  increases  with  increasing  Ohd. 

Other  researchers  have  examined  the  evolution  of  liquid  drops  in  another  liquid  with  a  density 
comparable  to  the  drop  density  moving  due  to  gravity.  The  deformation  of  miscible  liquid  drops 
at  low  Reynolds  numbers  was  studied  by  Kojima,  Hinch,  and  Acrivos30,  who  observed  that  the 
drops  form  vortex  rings.  The  stability  of  drops  moving  in  immiscible  fluids  was  investigated  by 
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Koh  and  Leal31’32,  who  showed  computational  results  for  zero  Reynolds  number  using  a  boundary 
integral  method  and  experimental  results  for  low  Reynolds  numbers.  A  similar  investigation  of  the 
instability  of  drops,  also  using  a  boundary  integral  method  was  reported  by  Pozrikidis33.  Experi¬ 
ments  by  Baumann,  Joseph,  Mohr,  and  Renardy34  showed  that  vortex  rings  can  also  be  created  in 
immiscible  liquids. 

A  few  investigators  have  simulated  the  deformation  and  breakup  of  liquid  drops  numerically. 
However,  due  to  the  difficulties  in  dealing  with  large  deformation  of  the  interface  and  in  accurately 
including  surface  tension  along  with  viscous  and  inertial  forces,  such  numerical  simulations  have 
often  been  based  on  considerable  simplifications.  The  steady  motion  of  deformable  axisymmetric 
drops  was  investigated  by  Dandy  and  Leal35  at  several  Reynolds  and  Weber  numbers  using  a 
finite-difference  method.  The  steady  rise  of  an  axisymmetric  drop  in  an  unbounded  surrounding 
fluid  was  examined  by  Volkov36  for  intermediate  Reynolds  numbers.  Bozzi,  Feng,  Scott,  and 
Pearlstein37  presented  finite  element  simulations  of  the  steady  motion  of  axisymmetric  drops  in 
bounded  domains.  Fritts,  Fyre,  and  Oran38  applied  a  two-dimensional  Lagrangian  finite  difference 
method  to  simulate  the  breakup  of  fuel  droplets  and  Liang,  Eastes,  and  Gharakhari39  presented 
simulations  of  axisymmetric  drop  breakup  using  a  Volume-of-Fluid  method  for  a  limited  number 
of  cases.  Other  numerical  studies  of  the  deformation  and  breakup  of  two-dimensional  drops  can 
be  found  in  Deng  and  Jeng40,  Deng,  Liaw,  and  Chou41,  Seung,  Chen,  and  Chen42,  and  Zaleski,  Li, 
and  Sued43.  However,  these  numerical  results  are  still  preliminary. 

In  spite  of  the  progress  made  by  previous  investigators,  several  aspects  of  the  secondary 
breakup  are  still  not  well  understood,  including  the  breakup  of  drops  at  high  pressure  and  tem¬ 
perature,  where  experimental  difficulties  are  encountered.  It  is  also  necessary  to  more  closely 
examine  the  time  dependent  characteristics  of  the  breakup.  In  existing  spray  models,  the  drop 
breakup  is  considered  to  occur  instantaneously.  Recent  experimental  evidence  indicates,  however, 
that  secondary  breakup  occurs  over  a  finite  amount  of  time.  Therefore,  it  is  possible  that  the  drop 
breakup  should  be  treated  as  a  time-dependent  process. 

In  this  paper,  a  numerical  method  based  on  a  front  tracking  technique  that  can  accommodate 
large  deformation  of  the  drops  is  developed  to  simulate  the  breakup  of  liquid  drops  accelerated  by  a 
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constant  body  force.  The  governing  equations  for  axisymmetric  geometry  are  solved  numerically 
on  a  non-uniform  grid  using  a  finite  difference  method.  The  drop  interface  is  represented  by 
connected  marker  points,  whose  positions  are  updated  explicitly  at  each  time-step.  Results  for  a 
wide  range  of  parameters  are  presented,  and  the  physical  significance  of  the  results  is  discussed. 


H.  FORMULATION  AND  NUMERICAL  METHOD 


The  physical  problem  and  computational  domain  are  sketched  in  Figure  1;  the  left  boundary 
is  the  axis  of  symmetry.  We  follow  the  motion  of  the  fluids  both  inside  and  outside  the  drop 
and  write  a  single  set  of  equations  for  the  whole  flow  field,  using  the  conservative  form  of  the 
governing  equations  to  allow  the  density  and  viscosity  to  change  discontinuously.  Surface  tension 
is  added  as  a  delta  function  to  provide  the  proper  interface  boundary  conditions.  Written  for  an 
axisymmetric  coordinate  system,  the  Navier-Stokes  equations  are: 


dpv  1  drpuv  dpv2  _ 

dt  **"  r  dr  dz 
dp  Id 

~Tz + 


dv  du 
dr  dz 


—  /  <7K,n5(x  —  Xf)dS  •  iz  —  paz 
Js 


Here,  u  and  v  are  the  velocity  components  in  the  radial  and  axial  directions,  p  is  the  pressure,  and 
p  and  p  are  the  discontinuous  density  and  viscosity  fields,  a  is  the  surface  tension,  k  is  twice  the 
mean  curvature,  ir  and  i2  are  the  radial  and  axial  components  of  the  surface  unit  normal  vector 
pointing  outward  from  the  drop,  and  5  is  a  three-dimensional  delta  function.  In  (1)  and  (2),  the 
surface  tension  is  treated  as  a  body  force.  The  integral  over  the  surface  of  the  drop,  S,  results  in 
a  force  that  is  smooth  and  continuous  along  the  drop  surface.  In  the  numerical  method,  the  delta 
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function,  8,  is  approximated  by  a  smooth  function  with  a  compact  but  finite  support.  The  constant 
acceleration  gives  rise  to  a  body  force  in  the  axial  direction  denoted  by  paz. 

The  above  equations  are  supplemented  by  the  incompressibility  condition: 

1  dru  dv 

“5-  +  =  0 

r  or  oz 

which,  when  combined  with  the  momentum  equations,  leads  to  an  elliptic  equation  for  the  pres- 


V  ■  j  =  -R 


where  R  is  the  divergence  of  the  vector  form  of  the  momentum  equations  (1)  and  (2),  excluding 
the  pressure  term. 

We  also  have  equations  of  state  for  the  physical  properties  of  the  drop  and  the  surrounding 


A  =  0  ;  A  =  0  (5) 

Dt  Dt 

where  D/Dt  is  the  material  derivative.  These  two  equations  state  that  the  physical  properties  of 
each  fluid  remain  constant. 

Dimensional  analysis  shows  that  four  independent  dimensionless  parameters  govern  the  dy¬ 
namics  of  drop  deformation  and  breakup.  When  the  drop  is  subject  to  an  acceleration  by  a  con¬ 
stant  body  force,  it  is  convenient  to  use  the  Eotvos  number,  Eo  (interchangeably  called  the  Bond 
number,  Bo)  and  the  Ohnesorge  number  of  the  drop,  Ohd,  defined  as: 


Eo  = 


azApD2 


Ohd  = 


Pd 

\f~pdPo 


where  A p  is  the  density  difference  between  the  drop  and  the  surrounding  fluid  and  D  is  the  initial 
diameter  of  the  drop.  The  density  and  the  viscosity  ratios: 

Pd  .  pd 
Po  ’  Po 
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can  be  selected  as  the  other  two  parameters.  The  viscosity  ratio  is  sometimes  replaced  by  the 
Ohnesorge  number  based  on  the  properties  of  the  surrounding  fluid: 


The  subscripts,  d  and  o,  denote  the  properties  of  the  drop  and  the  surrounding  fluid,  respectively. 
Time  is  non-dimensionalized  with  respect  to  the  drop  diameter  and  the  acceleration: 


For  drops  subject  to  impulsive  acceleration,  the  Weber  number,  We  defined  by: 


We  = 


p0U2D 


(10) 


(ID 


is  used  in  place  of  the  Eotvos  number  and  the  Ohnesorge  number  is  often  replaced  by  the  Reynolds 
number  defined  by: 


Re  = 


PoUD 

p0 


(12) 


Here,  U  is  the  initial  relative  velocity  between  the  drop  and  the  surrounding  fluid.  Time  is  non- 
dimensionalized  with  the  diameter  and  the  initial  relative  velocity: 


The  numerical  technique  used  for  the  simulations  presented  here  is  based  on  the  front¬ 
tracking/finite  difference  method  discussed  in  Unverdi  and  Tryggvason44.  The  code  employed 
in  the  present  study  is  an  axisymmetric  version  of  the  method.  Since  the  axisymmetric  code  runs 
much  faster  than  the  fully  three-dimensional  version,  it  allows  more  runs  and  higher  resolution. 
To  improve  the  efficiency  of  the  computations,  the  method  was  implemented  on  stretched  grids  to 
allow  clustering  of  grid  points  in  specific  regions. 

The  momentum  equations  and  the  continuity  equation  are  discretized  using  an  explicit  second- 
order  predictor-corrector  time-integration  method  and  a  second-order  centered  difference  approx¬ 
imation  for  the  spatial  derivatives.  The  discretized  equations  are  solved  on  a  fixed,  staggered  grid 
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using  the  Marker-and-Cell  method  developed  by  Harlow  and  Welch45.  The  full-slip  boundary 
condition  is  applied  to  all  four  boundaries. 

To  maintain  a  well  defined  boundary  between  the  drops  and  the  surrounding  fluid,  the  boundary 
is  marked  by  connected  points  (the  front)  that  are  advected  by  the  fluid  velocity,  interpolated  from 
the  fixed  grid.  The  new  position  of  the  marker  points  is  used  to  construct  a  new  density  field 
by  distributing  the  density  jump  to  the  grid  points  next  to  the  front  using  area  weighting,  and 
integrating  the  jump  to  find  the  density  everywhere.  Once  the  density  is  known,  the  viscosity  is 
set  as  a  function  of  the  density.  The  marker  points  are  also  used  to  find  the  surface  tension,  which 
is  then  assigned  to  the  nearest  grid  points  in  the  same  way  as  the  density  jump,  and  added  to  the 
discrete  Navier-Stokes  equations.  For  a  more  detailed  description  of  the  front  tracking  method, 
see  Unverdi  and  Tryggvason44  and  Tryggvason,  Bunner,  Ebrat,  and  Tauber46. 

The  implementation  of  the  numerical  technique  to  the  drop  breakup  problem  is  straightforward 
and  the  method  works  well  for  a  broad  range  of  parameters.  However,  as  pdj p0  is  raised,  the  com¬ 
putational  cost  increases,  partly  because  of  the  appearance  of  the  coefficient  1/p  in  the  pressure 
equation  (4)  but  also  because  the  effect  of  the  surrounding  fluid  is  weaker  when  the  density  ratio 
is  large  and  the  drop  travels  a  longer  distance  before  breaking  up.  In  order  to  avoid  having  to 
use  a  very  long  computational  domain,  we  move  the  computational  domain  with  the  drop.  The 
motion  of  the  domain  is  determined  from  the  solution,  and  an  extra  acceleration  term  is  added  to 
the  governing  equations  to  account  for  the  time  dependent  motion  of  the  domain.  The  boundary 
conditions  have  also  been  modified  to  include  a  constant  inflow  at  the  bottom  and  a  zero  velocity 
gradient  in  the  normal  direction  at  the  top. 

The  majority  of  the  simulations  presented  here  were  earned  out  on  HP  9000  workstations.  A 
typical  run  required  between  4000  and  120000  timesteps  and  took  12-240  hours,  depending  on 
the  parameters  of  the  problem. 

To  address  to  what  extend  the  observed  drop  evolution  can  be  described  by  an  inviscid  model, 
a  few  simulations  were  done  using  a  vortex  method.  The  interface  separating  the  drop  and  the 
surrounding  fluid  is  a  vortex  sheet  and  an  evolution  equation  for  the  vortex  sheet  strength  7  = 
(ud  —  u0)  •  s  can  be  derived  by  subtracting  the  tangential  components  of  the  Euler’s  equations  on 


8 


either  side  of  the  interface.  The  resulting  equation  is  (see,  for  example,  Tryggvason47): 

Dj  dV  ntfd\5  ,  iarl\  °  dK  nA\ 

_2+7_.s  =  2^(_.s  +  §_)--t--  (14) 

Here,  A  =  (pd  +  Po)/pd  +  Po  is  the  Atwood  number,  U  =  (1/2)  (ud  +  u0)  is  the  average  of  the 
velocities  on  either  side  of  the  vortex  sheet,  and  k  is  the  mean  curvature. 

Given  the  vortex  sheet  strength  7,  the  velocity  is  found  by  the  Biot-Savart  integral.  For  com¬ 
putational  purpose,  the  axisymmetric  vortex  sheet  is  discretized  by  a  finite  number  of  vortex  rings. 
The  azimuthal  integration  can  be  done  analytically  and  the  integral  is  therefore  replaced  by  a  sum¬ 
mation  over  the  discrete  vortex  rings.  The  radial  and  axial  velocity  at  point  on  ring  j  is  given 
by 

*  ■ ;  § r'  **  + 2K{k ,j)l  (15) 


*  -  +r-V~)  ^ + 2r,K{k,i)]  <16) 

Here  K(k)  is  the  Complete  Elliptic  integral  of  the  first  kind,  E(k)  is  the  Complete  Elliptic  integral 
of  the  second  kind,  and 

ktj  =  .  ■-  (17) 

\J irl  +  Tj)2  +  (Z3  -  Zl) 

The  Elliptic  integrals  can  be  computed  efficiently  by  a  polynomial  approximation.  When  the 
axisymmetric  vortex  sheet  is  replaced  by  discrete  vortex  rings,  the  rings  must  be  given  a  finite  core 
size  to  avoid  infinite  self-induced  velocity.  This  can  be  accomplished  simply  by  replacing  kij  by 

k,  =  (18) 

x/(r,  +  r,)2  +  (zJ -ziP+'S2 

where  6  is  a  small  regularization  parameter.  In  the  limit  of  iV  ->  00  and  5  — *  0  the  solution  will 
be  independent  of  the  exact  value  of  S  (except  at  isolated  points  where  roll-up  takes  place). 


III.  ACCELERATION  BY  A  CONSTANT  BODY  FORCE 


Numerical  simulations  are  presented  first  for  pd/ p0  =  10,  using  a  moving  computational  do¬ 
main.  To  examine  the  effect  of  the  density  ratio,  simulations  are  also  carried  out  for  a  small  density 


ratio,  Pd./ Po  =  1.15,  using  a  fixed  computational  domain.  For  each  density  ratio,  the  effects  of 
varying  the  other  dimensionless  parameters,  Eo,  Oha,  and  Ohd,  are  studied. 

A.  Validation 

In  order  to  validate  the  numerical  method,  grid  refinement  tests  were  performed.  Typical  results 
are  presented  in  Figure  2  where  the  shape  of  the  drop  is  plotted  at  time  intervals  At*  =  3.873,  using 
two  different  grids:  256  x  512  (left)  and  512  x  1024  (right).  The  non-dimensional  parameters 
are  pd/p0  -  1.15,  Eo  =  144,  Oh0  =  0.05,  and  Ohd  =  0.0466.  Initially  (f*  =  0),  the  drops 
are  spherical  and  the  velocities  are  zero  everywhere.  Despite  the  large  deformation  of  the  drop, 
the  results  agree  well.  In  Figure  3,  the  aspect  ratio  and  the  centroid  velocity  are  plotted  versus 
non-dimensional  time.  The  aspect  ratio  is  defined  as  the  maximum  width  of  the  drop  divided  by 
its  thickness  at  the  centerline.  The  centroid  velocity  is  found  by  taking  the  volume  average  of 
the  vertical  velocity  inside  the  drop.  The  results  corresponding  to  Figure  2  are  shown  along  with 
results  using  a  coarser  grid:  128  x  256.  The  result  from  the  128  x  256  grid  shows  a  small  difference 
but  the  two  finer  grids  give  nearly  identical  results. 

In  addition,  we  have  compared  our  results  to  the  steady  state  results  for  a  single  axisymmetric 
deformable  drop  computed  by  Dandy  and  Leal35.  They  specified  the  Reynolds  number  and  the 
Weber  number  and  found  the  drag  coefficient,  Cd,  as  a  part  of  the  solution.  In  our  transient 
simulation,  it  is  not  possible  to  specify  Re  and  We  a  priori,  since  the  velocity  of  the  drop  is 
computed  as  part  of  the  solution.  However,  once  the  drag  coefficient  is  known,  the  Eotvos  number 
and  Ohnesorge  number  can  be  found  by: 

Eo  =  hweCd ;  Oh0  =  (19) 

For  a  drop  translating  at  a  Reynolds  number  equal  to  100  and  Weber  number  equal  to  4,  with 
Pd/Po  =  0-91  and  Pd/Po  =  1,  Dandy  and  Leal35  found  Cd  =  0.919  in  an  unbounded  domain.  This 
gives  Eo  =  2.75  and  Oh0  =  0.02.  Our  computation  was  done  using  a  256  x  768  grid  in  a  domain 
5  and  15  times  the  initial  diameter  of  the  drop  in  the  radial  and  axial  directions,  respectively. 
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The  computed  Reynolds  and  Weber  numbers  differ  from  those  given  by  Dandy  and  Leal35  by 
less  than  1%  when  the  drop  reaches  steady  state.  For  Pd/Po  =  4  and  the  same  Re,  We,  pd/p0 
as  in  the  previous  case,  Dandy  and  Leal35  found  Cd  =  1.10.  Our  computation  was  done  using 
Eo  =  3.3,  Oh0  =  0.02,  and  the  same  pd/p0  and  pd/p0,  with  the  same  resolution  and  domain  size 
as  in  the  previous  case.  The  result  gives  Cd  =  1.13,  which  is  approximately  3%  higher.  This  is 
due  to  the  finite  size  of  our  computational  domain.  Computations  using  domains  of  half  and  twice 
the  original  size  in  the  radial  direction  yield  Cd  =  1-19  and  1.12,  respectively. 

B.  The  Boussinesq  approximation 

Before  presenting  further  computational  results,  we  pause  to  examine  the  validity  of  the 
Boussinesq  approximation.  The  Boussinesq  approximation  states  that  if  the  density  difference 
is  small,  density  variations  are  only  important  when  multiplied  by  gravity.  A p  is  therefore  no 
longer  an  independent  parameter  and  it  is  sufficient  to  simulate  the  breakup  for  only  one  value  of 
the  density  ratio  in  this  limit.  Results  for  other  values  of  A p  can  be  obtained  by  simply  rescaling 
time.  For  a  discussion  of  the  Boussinesq  approximation  to  stratified  flows,  see,  for  example,  Dahm, 
Scheil,  and  Tryggvason48.  The  relative  magnitude  of  the  density  difference  is  better  expressed  by 


the  Atwood  number,  defined  by: 


P  d  ~  Po 
Pd  +  Po 


When  A  is  sufficiently  small  (pd/p0  ->  1),  time  and  velocities  can  be  scaled  by  the  average  static 
pressure  to  yield: 


DKAa* 


CAaJ)' 


1  Aa.D 


The  Eotovs  number  and  the  Ohnesorge  number  must  also  be  redefined  as: 


Eo  — 


PavAdzD 
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(24) 


1^0 

y/ PavDa 


where  pav  =  0.5 {pd  +  p0).  Note  that  the  constant  acceleration  appears  only  as  Aaz  instead  of  az 
alone. 


In  order  to  check  the  validity  of  the  Boussinesq  approximation,  tests  were  done  for  different 
values  of  A.  In  Figure  4,  results  for  Eo  =  72,  Oh0  =  0.241,  and  pd/Po  =  1  are  presented. 
The  computation  was  done  using  a  256  x  768  grid  and  the  size  of  the  computational  domain  was 
5  and  15  times  the  initial  drop  diameter  in  the  radial  and  the  axial  direction,  respectively.  The 
aspect  ratio  a  and  the  non-dimensionalized  centroid  velocity  Vc  are  plotted  versus  t  in  (a)  and 
(b),  respectively.  In  each  graph,  simulations  using  four  values  of  A  are  shown:  0.07,  0.11,  0.2, 
and  0.33,  corresponding  to  the  density  ratios:  1.15, 1.25, 1.5,  and  2.0.  The  plots  confirm  that  the 
scaling  works  well  when  A  is  less  than  about  0.2  to  0.3.  The  ability  to  cover  this  density  range  by 
a  single  simulation  is  obviously  a  considerable  simplification. 


C.  Effect  of  Eo  at  small  Oh 


When  Oh  is  small  and  surface  tension  is  much  more  important  than  viscous  stresses.  Oh  has 
little  influence  on  the  breakup  and  Eo  is  the  only  controlling  parameter.  Here,  we  present  results 
for  different  Eo  when  Oh  is  small.  When  a  drop  is  set  into  motion  by  a  constant  body  force,  the 
hydrodynamic  pressure  is  higher  at  the  poles  and  lower  at  the  equator  and  the  drop  deforms  into 
an  oblate  ellipsoid.  This  deformation  is  opposed  by  the  surface  tension.  Depending  on  the  relative 
strength  of  the  pressure  forces  and  the  surface  tension,  measured  by  Eo,  different  breakup  modes 
are  observed. 

In  Figure  5,  the  effect  of  Eo  is  presented  for  Pd/po  —  10  ar*d  Oh0  =  Ohd  —  0.05.  The 
simulations  are  done  using  a  moving  coordinate  system  where  the  origin  is  fixed  at  the  centroid 
of  the  drop.  The  domain  has  dimensions  of  five  and  fifteen  times  the  initial  drop  diameter  in  the 
radial  and  axial  directions,  respectively.  The  centroid  of  the  drop  is  fixed  at  a  position  five  times 
the  initial  drop  diameter  above  the  bottom  boundary.  The  evolution  of  the  drop  is  shown  for  nine 
different  values  of  Eo  (a  to  i).  In  each  column,  the  drop  interface  is  plotted  at  fixed  time  intervals. 
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The  separation  between  two  successive  drops  is  equal  to  the  distance  that  the  drop  travels  during 
the  time  interval. 

In  (a),  the  drop  is  shown  for  Eo  =  12.  As  the  drop  starts  falling,  the  back  side  becomes  flat 
while  the  front  side  retains  a  rounded  shape.  After  the  initial  deformation,  the  drop  reaches  a  steady 
state  and  no  further  change  in  the  drop  shape  is  seen.  When  Eo  is  increased  to  24  in  (b),  the  drop 
deformation  is  more  pronounced.  Initially,  the  drop  assumes  a  shape  similar  to  that  shown  in  (a), 
but  then  the  back  of  the  drop  becomes  increasingly  more  convex  and  eventually  the  drop  deforms 
into  a  thin  disk-like  shape  that  moves  at  a  nearly  steady  state.  The  drop  shown  in  (c)  for  Eo  =  28.8 
evolves  in  the  same  way  until  it  has  deformed  into  a  disk-like  shape.  Then  the  thickness  of  the 
drop  near  the  symmetry  axis  continues  to  decrease,  and  most  of  the  drop  fluid  moves  outward 
toward  the  edge  of  the  drop.  Finally,  the  center  of  the  front  surface  is  pushed  upward,  forming  a 
backward-facing  bag.  At  this  stage,  most  of  the  drop  fluid  is  contained  in  the  annular-shaped  rim. 
As  time  progresses,  the  bag  expands  both  radially  outward  and  vertically  upward.  Experimental 
evidence  indicates  that  the  drop  will  eventually  break  into  small  drops.  The  evolution  shown  in  (d) 
for  Eo  =  36  is  very  similar  to  that  in  (c),  displaying  a  backward-facing  bag.  The  only  difference 
is  that  the  rate  of  deformation  is  higher  and  the  backward-facing  bag  expands  more  rapidly. 

When  Eo  is  further  increased  to  48  in  (e),  a  different  mode  of  breakup  is  observed.  The  initial 
deformation  is  not  very  different  from  the  previous  cases,  and  an  indentation  develops  on  the  back 
surface,  but  instead  of  deforming  into  a  disk-like  shape,  the  drop  remains  relatively  thick  near  the 
symmetry  axis  and  the  edge  of  the  drop  is  swept  back  in  the  downstream  direction.  A  large  wave 
then  develops  on  the  drop  interface  and  as  this  wave  propagates,  the  drop  deforms  in  an  erratic 
manner.  The  evolution  of  the  drop  shown  in  (f)  for  Eo  =  60  reveals  another  mode  of  deformation. 
The  initial  evolution  is  similar  to  the  previous  cases,  but  the  results  are  different  at  later  times.  As 
the  indentation  at  the  top  progressively  deepens,  the  drop  does  not  deform  into  a  thin  disk-like 
shape.  Instead,  the  edge  of  the  drop  is  deflected  in  the  downstream  direction  and  drawn  out  into  a 
thin  film  with  a  blob  of  drop  fluid  at  the  end.  The  appearance  of  this  film  is  similar  to  the  skirted 
drop  shapes  observed  in  experimental  studies  of  liquid  drops  moving  at  steady  state  (Wairegi  and 
Grace49).  The  center  portion  of  the  drop,  however,  maintains  a  convex  shape  and  its  thickness  at 
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the  symmetry  axis  stops  decreasing.  Similar  drop  deformations  are  observed  at  even  higher  Eo- 
72, 96,  and  144  as  shown  in  (g),  (h),  and  (i),  respectively. 

Based  on  these  results,  the  evolution  of  drops  with  pd/ p0  =  10  at  a  small  Oh  can  be  classified 
into  four  categories  in  order  of  increasing  Eo:  steady  deformation,  formation  of  a  backward-facing 
bag,  transient  breakup  with  a  complex  shape,  and  stripping  or  shearing  of  a  film  from  the  edge  of 
the  drop.  It  is  evident  from  Figure  5  that  drops  breaking  up  in  the  backward-facing  mode  travel 
a  much  longer  distance  than  those  breaking  up  in  the  shear  breakup  mode.  Also  note  that  for  the 
same  breakup  mode,  the  rate  of  drop  deformation  increases  as  Eo  increases. 

In  Figure  6,  the  evolution  of  a  drop  with  a  small  density  ratio,  Pd/po  =  1.15,  is  shown  for 
different  Eo.  Again,  values  of  Ohnesorge  numbers,  Oh0  =  0.05  and  Ohd  =  0.0466  are  chosen  so 
that  viscous  stresses  are  small  compared  to  surface  tension.  The  computations  were  done  using  a 
fixed  coordinate  system.  When  Eo  is  small,  the  drop  deforms  into  an  oblate  ellipsoid  and  moves 
with  a  steady  state  shape  as  shown  in  (a)  for  Eo  =  12.  When  Eo  is  increased  to  24  (b),  the  drop 
deforms  more  and  eventually  forms  a  backward-facing  bag  as  observed  for  pd/ p0  =  10  in  Fig¬ 
ure  5(c).  In  (c),  where  Eo  is  48,  the  drop  moves  with  an  essentially  steady  convex  shape,  showing 
no  sign  of  bag  formation.  Compared  to  its  high  density  ratio  counterpart  shown  in  Figure  5(e),  the 
overall  deformation  is  reduced.  When  Eo  is  96  in  (d),  the  indentation  at  the  back  of  the  drop  deep¬ 
ens  continuously  until  it  reaches  the  front  of  the  drop,  creating  a  forward-facing  bag.  Eventually, 
however,  the  heavier  edge  falls  faster  than  the  thin  bag.  This  formation  of  a  forward-facing  bag  is 
different  from  the  shear  breakup  mode  observed  in  Figure  5(f)-(i),  where  a  significant  portion  of 
the  fluid  remains  near  the  symmetry  axis  while  a  thin  film  is  pulled  away  from  the  edge.  In  (e), 
Eo  is  further  increased  to  144.  The  overall  evolution  is  similar  to  (d),  but  the  rate  of  deformation 
is  slightly  faster. 

In  Figure  7,  vorticity  contours  (left)  and  streamlines  with  respect  to  a  frame  moving  with 
the  drop  (right)  are  plotted  at  a  few  selected  times  for  the  drop  shown  in  Figure  5(c).  Most  of 
the  vorticity  is  created  at  the  outer  edge  of  the  drop,  as  expected,  and  the  streamlines  show  the 
formation  of  a  large  wake  behind  the  drop.  The  pressure  difference  between  the  front  stagnation 
point  and  the  wake  causes  the  formation  of  the  backward-facing  bag. 
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Figure  8  shows  vorticity  contours  (left)  and  streamlines  (right)  at  a  few  selected  times  for  the 
drop  shown  in  Figure  5(f).  Although  vorticity  generated  at  the  drop  surface  accumulates  into  a 
large  wake  as  in  Figure  7,  the  more  deformable  drop  is  continuously  deformed  by  the  flow  (as 
seen  by  streamlines  crossing  its  boundary)  and  the  edge  is  pulled  back  by  the  flow.  In  the  last 
frame,  the  large  wake  formed  initially  separates  from  the  drop,  leaving  the  film  to  move  toward 
the  axis  due  to  the  flow  around  the  smaller  remaining  drop. 

In  Figure  9,  vorticity  contours  (left)  and  streamlines  (right)  are  shown  for  the  drop  shown  in 
Figure  6(e).  Here,  the  vorticity  generated  at  the  interface  moves  with  the  drop,  forming  a  dipole 
that  continuously  deforms  the  interface  into  a  forward  facing  bag. 

Figures  10  and  11  show  the  centroid  velocity  of  the  drop  Vc  plotted  versus  t*  for  the  drops 
shown  in  Figures  5  and  6,  respectively.  Since  the  velocity  is  non-dimensionalized  by  y/a^D,  the 
graphs  for  different  values  of  Eo  all  have  the  same  initial  slope.  After  the  initial  acceleration,  the 
drop  deformation,  which  depends  on  Eo,  determines  the  velocity.  In  Figure  10,  drops  with  low  Eo 
deform  less  and  therefore  move  faster  than  drops  with  high  Eo.  The  lowest  Eo  drop  ( Eo  =  12) 
asymptotically  reaches  a  steady  state  velocity,  but  the  other  drops  all  slow  down  as  they  start 
deforming.  The  Eo  =  24  drop  also  reaches  a  steady  velocity.  The  drops  that  undergo  bag  breakup 
first  behave  like  the  Eo  =  24  drop,  but  as  the  bag  forms,  the  drops  slow  down  rapidly.  At  the  very 
end,  the  Eo  =  36  drop  speeds  up  again,  as  the  rim  of  the  drop  starts  falling  independently  of  the 
bag.  The  rest  of  the  drops  all  slow  down  rapidly  as  they  are  stretched  perpendicular  to  the  flow,  and 
all  speed  up  again  as  the  thin  film  pulled  from  their  edges  folds  back  toward  the  axis.  The  results 
for  the  small  density  difference  in  Figure  11  show  a  similar  trend,  but  with  a  few  differences.  The 
transient  drop  (Eo  =  48)  reaches  a  velocity  that  is  nearly  the  same  as  the  velocity  of  the  drop 
moving  with  a  steady  deformed  shape  (Eo  =12)  and  the  reduction  in  speed  due  to  bag  breakup  is 
smaller  than  in  Figure  10. 

In  Figures  12  and  13,  the  surface  area  S  (normalized  by  the  initial  value  S0)  is  plotted  versus  t* 
for  the  drops  shown  in  Figures  5  and  6,  respectively.  The  graphs  for  pdj p0  =  10  in  Figure  12,  show 
that  a  rapid  increase  of  the  surface  area  takes  place  when  breakup  occurs,  and  that  the  backward¬ 
facing  bag  breakup  mode  takes  longer  than  the  shear  breakup  mode.  The  drops  undergoing  a  shear 
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breakup  show  a  reduction  in  surface  area  when  the  film  moves  toward  the  symmetry  axis.  The 
graphs  for  pd/p0  -  1.15  in  Figure  13  also  display  a  rapid  increase  of  the  surface  area  when  the 
drops  break  up.  However,  the  drops  with  the  highest  Eo  show  a  rapid  increase  in  surface  area  after 
the  rim  starts  falling,  and  the  surface  area  of  the  drop  undergoing  bag  breakup  grows  relatively 
slowly  compared  to  the  higher  density  ratio  drops. 

D.  Effect  of  Oh 

Figure  14  illustrates  the  effect  of  the  Ohnesorge  number  (Oh-the  non-dimensional  viscosity) 
for  drops  with  a  finite  density  ratio,  pd/pa  =  10.  The  drops  are  shown  at  several  stages.  Here, 
Ohd  is  equal  to  Oha.  The  case  where  Ohd  is  different  from  Oh0  will  be  discussed  in  the  next 
section.  In  the  top  row  (a-c),  Oh  =  0.05,  0.125,  and  0.25,  from  left  to  right,  and  Eo  =  28.8.  The 
Oh  =  0.05  case  (a)  has  already  been  shown  in  Figure  5(c),  but  is  included  here  for  comparison. 
The  initial  deformation  of  all  three  drops  is  similar,  but  whereas  the  Oh  =  0.05  drop  (a)  deforms 
into  a  backward  facing  bag,  the  other  two  drops  reach  a  steady  state  shape.  Of  those,  the  less 
viscous  drop  (b)  is  flatter. 

In  the  bottom  row  (d— f),  Eo  is  increased  to  144  and  the  evolution  of  the  drops  is  presented 
for  the  same  three  values  of  Oh  as  in  the  top  row.  In  (d),  the  drop  already  shown  in  Figure  5(i)  is 
included  for  reference.  This  drop  undergoes  a  so-called  shear  (or  boundary  stripping)  breakup.  The 
Oh  =  0.125  drop  (e)  shows  a  similar  evolution  as  the  drop  in  (d),  although  the  rate  of  deformation 
is  reduced  slightly.  The  center  portion  of  the  drop  still  contains  a  significant  amount  of  drop  fluid 
and  formation  of  a  backward-facing  bag,  which  requires  the  formation  of  a  very  thin  film  of  fluid 
near  the  symmetry  axis,  does  not  occur.  In  contrast,  the  center  portion  of  the  drop  in  (f)  is  drained 
completely,  and  the  drop  forms  a  backward  facing  bag. 

Based  on  the  results  shown  in  Figure  14,  it  is  clear  that  increasing  both  Oh0  and  Ohd  simulta¬ 
neously  results  in  a  translation  of  the  boundaries  between  the  breakup  modes  to  higher  Eo. 

Figure  15  illustrates  the  effect  of  viscosity  on  the  initial  deformation  of  drops  with  pd/ p0  = 
1.15.  In  addition  to  runs  with  a  finite  viscosity,  we  show  simulations  with  zero  viscosity,  obtained 
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by  an  axisymmetric  vortex  method  (see  Dahm,  Frieler,  and  Tryggvason50).  In  the  top  row,  Eo  —  24 
and  Oh0  is  0.05,  0.025,  0.01  and  0  from  left  to  right.  In  all  cases,  pd/Ho  =  1.  While  the  initial 
acceleration  is  dominant,  all  the  drops  evolve  in  the  same  way.  As  time  progresses,  viscosity 
effects  become  important  and  the  viscous  drops  will  eventually  develop  a  backward-facing  bag 
due  to  the  formation  of  a  wake.  See  Figure  6(b)  for  further  deformation  of  the  drop  in  (a).  The 
inviscid  drop  (d),  on  the  other  hand,  loses  fluid  to  a  film  pulled  off  its  edge.  Since  this  drop  does 
not  develop  a  vortical  wake,  it  does  not  form  a  backward  facing  bag.  Similar  deformation  is  also 
seen  for  inviscid  bubbles51 . 

In  the  bottom  row,  the  limit  of  zero  surface  tension,  i.e.  Eo  —  oo  is  investigated.  Since  Oh0  is 
defined  with  the  surface  tension  in  the  denominator,  it  cannot  be  used  as  a  measure  of  viscosity  in 
this  case.  Instead,  another  non-dimensional  number  is  defined: 


where  Ar  is  the  Archimedes  number.  Results  for  three  values  of  p  =0.01021  (e),  0.00513  (f),  and 
0.00204  (g)  are  compared  with  results  for  ft  =  0  (h).  In  all  cases,  Hd/ Ho  =  1-  From  the  plots,  it 
can  be  seen  that  since  there  is  no  surface  tension  limiting  the  deformation,  all  the  drops  evolve  in 
a  similar  way:  first  an  indentation  forms  at  the  top,  then  the  drops  deform  into  a  forward-facing 
bag  with  a  thick  edge.  The  effect  of  fi  on  the  overall  shape  of  the  drop  is  relatively  small,  with  the 
exception  of  the  rollup  of  the  edge  which  increases  as  p,  is  reduced. 

E.  Effect  of  the  viscosity  ratio 

The  results  presented  so  far  are  all  for  drops  moving  in  another  fluid  that  has  the  same  or  almost 
the  same  viscosity.  The  effect  of  the  viscosity  ratio  is  shown  in  Figure  16,  where  the  drops  are 
shown  for  several  values  of  the  governing  parameter.  In  (a)  and  (b),  pd / p0  =  10  and  Oh0  =  0.05. 
Eo  is  72  in  (a)  and  144  in  (b).  In  each  row,  the  drop  shape  is  shown  at  a  fixed  t*  for  Ohd  =  0.05, 
0.125, 0.25,  and  1.25,  increasing  from  left  to  right.  The  evolution  of  the  drops  in  (a)  is  qualitatively 
similar  for  the  three  lower  values  of  Ohd.  They  all  show  a  shear  breakup  mode  in  which  a  film 
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of  drop  liquid  is  pulled  away  from  the  edge  of  the  drop.  The  drop  with  the  highest  Ohd  in  the 
rightmost  column,  however,  has  not  progressed  as  far,  and  will  eventually  form  a  forward-facing 
bag.  The  comparison  in  (b),  for  Eo  =  144,  shows  the  trend  observed  in  (a).  Comparisons  for 
drops  with  pd/p0  =  1.15  are  presented  in  (c)  and  (d).  In  (c),  Eo  =  24  and  Oh0  =  0.05.  The  drops 
are  shown  for  0hd=0.0093, 0.0466, 0.2331,  and  1.1656  (from  left  to  right)  at  t*  =  63.2.  The  drops 
with  the  three  lower  viscosity  ratios  form  a  backward-facing  bag  and  the  drop  deformation  is  most 
pronounced  when  the  viscosity  ratio  is  smaller.  In  contrast,  the  most  viscous  drop  develops  a 
steady  disk-like  shape.  In  (d),  Eo  =  144  and  Oh0  =  0.25  and  the  drops  are  shown  for  four 
different  values  of  the  drop  Ohnesorge  numbers,  Ohd  =  0.0093, 0.0466, 0.2331,  and  1.1656  (from 
left  to  right)  at  t*  =  27.1, 27.1, 46.5,  and  62.0,  respectively.  The  times  are  not  the  same  because 
as  the  viscosity  ratio  increases,  the  drops  deform  much  more  slowly.  Here,  the  edge  of  the  drop 
is  pulled  backward  into  a  thin  skirt  for  the  three  lower  viscosity  ratios.  The  most  viscous  drop, 
Ohd  =  1.1656,  does,  on  the  other  hand,  form  a  backward-facing  bag. 

In  Figure  17,  the  evolution  of  the  centroid  velocity  is  plotted  for  the  drops  shown  in  Figure  16. 
Initially,  while  the  drops  are  nearly  spherical,  acceleration  is  independent  of  the  viscosity  of  the 
drop  fluid.  As  the  drops  start  deforming,  they  slow  down  due  to  increased  drag.  For  the  drops  with 
Pd/po  =  10  shown  in  (a)  and  (b),  the  higher  viscosity  drops  deform  more  slowly  and  therefore 
move  faster.  At  a  later  time,  however,  the  most  viscous  drop  in  (a)  forms  a  backward  facing  bag 
and  slows  down  continually.  In  (c),  where  the  density  ratio  is  lower,  the  most  viscous  drop  reaches 
a  steady  state  shape  and  velocity.  The  other  drops  all  form  bags  and  are  relatively  unaffected  by 
changes  in  the  drop  viscosity.  In  (d),  the  low  viscosity  drops  speed  up  again,  once  a  skirt  has  been 
pulled  off  their  edges,  indicating  that  the  skirt  has  no  significant  effects  on  the  motion  at  this  stage. 
The  most  viscous  drop,  on  the  other  hand,  forms  a  backward  facing  bag  and  continues  to  slow 
down. 

Increasing  the  drop  viscosity  reduces  its  rate  of  deformation  and  in  some  cases,  this  can  result 
in  different  breakup  modes,  changing  a  shear  breakup  to  a  bag  breakup,  and  a  bag  breakup  to  a 
steady-state  shape. 
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F.  Deformation  and  breakup  regime  maps 


To  summarize  the  results  of  the  various  simulations,  deformation  and  breakup  maps  are  pre¬ 
sented  in  Figures  18  and  19.  In  the  maps,  we  mark  the  location  of  each  simulation  in  the  Ohd- 
Eo  plane,  using  a  different  symbol,  depending  on  the  deformation  and  breakup  mode  observed. 
Similar  breakup  maps  have  been  used  to  present  experimental  results  by  numerous  investigators. 
Figure  18  shows  the  result  of  simulations  with  Pd/Po  =  10..  Three  maps  are  shown  (a-c),  corre¬ 
sponding  to  different  ambient  Ohnesorge  numbers. 

The  map  for  a  relatively  small  Oh0  =  0.05  (a)  shows  that  increasing  the  magnitude  of  Eo 
at  a  fixed  drop  Ohnesorge  number  (Ohd  =  0.05)  results  in  the  following  transitions  between 
the  different  breakup  modes:  oblate  ellipsoid  ->•  backward-facing  bag  mode  transient  breakup 
shear  breakup  mode.  Changing  Ohd  for  a  fixed  Eo,  on  the  other  hand,  yields  only  minor 
differences  in  the  breakup  mode.  Increasing  Ohd  from  0.05  to  1.25,  when  Eo  is  fixed  at  28.8, 
for  example,  does  not  change  the  breakup  mode.  For  Eo  —  72  and  144,  a  change  from  a  shear 
breakup  mode  to  a  forward-facing  bag  mode  is  observed,  as  Ohd  is  increased  from  0.25  to  1.25. 
The  difference  between  these  two  breakup  modes  is,  however,  not  as  significant  as  that  between 
the  backward-facing  bag  and  the  shear  breakup  modes. 

When  Oh0  is  increased  to  0.125,  map  (b),  the  increased  viscosity  of  the  surrounding  fluid 
slows  the  drop  down  and  reduces  the  rate  of  deformation.  At  Ohd  =  0.05,  increasing  Eo  yields 
the  following  transitions  between  breakup  modes:  deformed  drop  — >  backward-facing  bag  — > 
shear  breakup.  When  the  boundaries  between  these  breakup  modes  are  compared  to  the  same 
Ohd  in  map  (a),  some  differences  are  observed.  The  backward-facing  bag,  which  was  observed 
at  Eo  =  28.8  in  (a),  is  now  seen  at  Eo  =  36.  At  Eo  =  48,  the  transient  breakup  is  no  longer 
observed  and  instead  we  see  a  backward-facing  bag. 

The  effect  of  changing  Ohd  at  a  fixed  Eo  is  also  examined  in  (b).  When  Eo  =  28.8,  changing 
Ohd  from  0.05  to  1.25  results  in  only  minor  differences.  The  Ohd  =  0.05  drop  displays  a  prolate 
shape  after  an  initial  oscillatory  motion  but  drops  with  higher  Ohd  deform  into  oblate  ellipsoids 
with  an  indentation  (or  a  dimple)  at  the  top.  At  Eo  =  72,  the  effect  of  changing  Ohd  is  more 
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significant.  Only  the  Ohd  =  0.05  drop  shows  a  shear  breakup  mode,  while  the  higher  Ohd  drops 
all  form  a  backward-facing  bag.  At  Eo  =  144,  on  the  other  hand,  the  breakup  modes  are  generally 
similar  to  those  observed  in  (a),  showing  no  significant  change  with  Ohd. 

When  Oh0  is  increased  to  0.25  in  map  (c),  the  effect  of  the  drop  viscosity  becomes  more 
significant.  At  Eo  =  28.8,  four  simulations  with  Ohd  ranging  from  0.05  to  1.25  show  deformed 
drops,  which  are  similar  to  those  observed  in  (b)  at  the  same  Eo.  At  Eo  =  72,  four  simulations 
with  Ohd  in  the  same  range  all  show  the  formation  of  a  backward-facing  bag.  This  is  different 
from  the  result  in  (b)  for  the  same  Eo,  where  the  shear  breakup  mode  changes  to  a  backward¬ 
facing  bag  mode  as  Ohd  increases.  At  Eo  =  144,  only  the  Ohd  =  0.05  drop  shows  the  shear 
breakup  mode  seen  in  (a)  and  (b).  The  Ohd  =  0.125  drop  deforms  into  a  forward-facing  bag  and 
the  higher  Ohd  drops  display  a  backward-facing  bag  mode. 

In  Figure  19,  deformation  and  breakup  regime  maps  are  presented  for  drops  with  pd/ p0  =  1.15. 
Three  maps  are  shown  in  (a)-(c),  for  ambient  fluids  with  Oh0  =  0.05, 0.25,  and  1.25.  The  map 
for  Oh0  =  0.05  (a),  displays  the  following  transitions  between  breakup  modes  as  Eo  is  increased: 
oblate  ellipsoid  -4  backward-facing  bag  -4  oscillating  indented  drop  -4  forward-facing  bag.  It  is 
clear  that  increasing  Ohd  has  no  major  effects.  The  only  exception  is  when  Ohd  becomes  large 
(>  1)  and  Eo  is  relatively  low.  In  this  case,  the  backward-facing  bag  breakup  mode  is  replaced  by 
a  steadily  moving  indented  drop. 

When  Oh0  is  increased  to  0.25  (b),  a  backward-facing  bag  mode  is  no  longer  observed  when 
Eo  =  24.  Instead,  a  steadily  moving  indented  drop  is  seen  for  the  Ohd  range  investigated.  The 
breakup  mode  at  Eo  =  144  also  changes  from  a  forward-facing  bag  mode  to  a  skirted  drop  when 
Ohd  is  small  (<  1).  A  forward-facing  bag  mode  is  observed  at  Eo  —  288.  As  in  (a),  no  noticeable 
effects  of  changing  Ohd,  at  a  fixed  Eo,  are  observed  as  long  as  Ohd  <  1.  When  Ohd  >  1,  the  drop 
develops  a  backward-facing  bag  when  Eo  =  144. 

In  map  (c),  Oh0  =  1.25  and  the  high  viscosity  prevents  nearly  all  deformation.  When  Eo  =  24, 
the  drop  remains  an  oblate  ellipsoid  but  for  Eo  —  144,  the  drops  develop  an  indentation  at  the  back. 
The  indentation  of  the  more  viscous  drop  (Ohd  =  1.1656)  deepens  continuously  until  it  reaches 
the  bottom  surface  of  the  drop,  forming  a  vortex  ring. 
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G.  Conclusions 


The  deformation  and  breakup  of  axisymmetric  drops,  accelerated  by  a  constant  body  force, 
have  been  studied  by  numerical  simulations.  Results  are  presented  for  two  density  ratios, 
pd/po  =  1.15  and  pd/p0  =  10.  For  the  lower  density  ratio,  the  Boussinesq  approximation  is  valid 
and  the  results  therefore  apply  for  other  low  density  ratios  by  the  simple  rescaling  discussed  in 
section  HI.  For  low  Ohnesorge  numbers,  the  Eotvos  number  and  the  density  ratio  are  the  main 
controlling  parameters.  At  low  density  ratios  the  drop  deforms,  but  does  not  break  up  for  Eo  less 
than  about  18.  For  18  <  Eo  <  36  (approximately),  the  drop  breaks  up  by  the  formation  of  a  back¬ 
ward  facing  bag.  Transient  breakup  is  observed  for  Eo  around  48,  and  at  Eo  larger  than  about  60, 
the  drop  evolves  into  a  forward  facing  bag. 

The  formation  of  a  forward-facing  bag  takes  place  very  quickly  (the  drop  has  moved  only  3—4 
times  its  initial  diameter  when  the  bag  is  formed)  and  is  essentially  an  inviscid  phenomenon.  The 
formation  of  a  backward-facing  bag,  on  the  other  hand,  takes  significantly  longer  (the  drop  has 
moved  8-10  times  its  initial  diameter).  A  comparison  with  results  obtained  by  an  inviscid  vortex 
method  shows  that  the  backward  facing  bag  is  a  viscous  phenomenon,  due  to  the  formation  of  a 
low  pressure  wake  behind  the  drop.  Furthermore,  the  surface  area  of  the  drop  increases  at  a  faster 
rate  in  the  forward-facing  bag  mode. 

As  Oh  is  increased,  the  effect  of  the  viscosity  reduces  the  rate  of  deformation.  At  low  Eo, 
while  the  drop  flattens,  its  center  does  not  drain  completely  and  backward-facing  bag  does  not 
form.  As  Eo  becomes  larger,  the  edges  of  the  drop  are  pulled  outward  and  sheared  off,  leading  to 
a  “skirted”  drop. 

When  Oh  becomes  very  large,  the  drop  deforms  into  an  oblate  ellipsoid  at  low  Eo.  At  high  Eo, 
baroclinically  generated  vorticity  causes  indentation  of  the  back  of  an  accelerated  drop,  but  rapid 
diffusion  of  vorticity  prevents  the  roll-up  observed  for  lower  Oh.  As  this  indentation  continues 
to  grow,  the  drop  finally  breaks  into  a  toms.  Similar  evolution  has  been  seen  in  simulations  of 
initially  oblate  drops  in  Stokes  flow  (Koh  and  Leal31,  Pozrikidis33). 

The  effect  of  the  viscosity  ratio  is  small  when  Oh0  is  small.  Although  there  are  differences  in 
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the  detailed  shape  of  the  drop,  the  overall  evolution  is  generally  similar  when  Ohd  is  varied  at  fixed 
Eo  and  OhB  unless  Ohd  is  very  large.  For  larger  Oh0,  the  effect  of  the  viscosity  ratio  becomes 
more  significant  and  as  Ohd  increases,  the  boundaries  between  the  different  breakup  modes  are 
moved  to  higher  Eo. 

At  higher  density  ratio  (pd/po  =  10)  the  evolution  is  similar  to  the  low  density  ratio  case  and 
the  effect  of  the  governing  parameters  is  also  similar.  There  are,  however,  a  few  important  differ¬ 
ences.  At  large  Eo,  the  forward  facing  bag  seen  for  the  low  density  case  is  replaced  by  a  shear 
breakup  mode  where  the  edge  of  the  drop  is  pulled  in  the  downstream  direction,  forming  a  blob  of 
drop  fluid  connected  to  the  main  drop  by  a  thin  film.  The  skirted  drop  in  the  low  density  ratio  case 
is  similar  in  shape  (except  that  no  blob  was  formed).  The  skirt,  however,  appears  only  at  a  rela¬ 
tively  higher  Oh  and  grows  slowly  once  it  has  formed.  In  contrast,  the  shear  breakup  mode  in  the 
higher  density  ratio  case  occurs  across  the  Oh  range  investigated  here.  The  boundaries  between 
the  breakup  modes  at  low  Oh  remain  essentially  the  same  as  for  the  low  density  ratio  case. 

In  most  practical  combustion  systems,  the  density  difference  between  the  liquid  fuel  and  the 
high  pressure  gas  is  considerably  smaller  than  at  atmospheric  temperature  and  pressure.  For  diesel 
engines,  pd/p0  =  32  -  53,  for  example,  (see  Heywood52)  and  pd/p0  of  order  unity  is  common 
in  rocket  motors.  Nearly  all  experimental  studies  of  secondary  breakup  of  drops,  however,  have 
been  done  at  atmospheric  pressures.  The  present  study  approaches  the  breakup  problem  from 
the  small  density  ratio  limits,  thus  complementing  previous  work.  Covering  the  gap  for  density 
ratios  between  those  studied  here  and  the  experiments  is  within  the  range  of  present  computational 
capabilities,  but  requires  considerably  longer  computational  times. 


IV.  IMPULSIVE  ACCELERATION 

The  breakup  of  drops  subject  to  impulsive  acceleration  has  been  examined  for  two  density 
ratios  (pd/p0  =  1.15  and  (pd/p0  =  10).  As  shown  in  the  last  section,  a  simple  rescaling  of  time 
allows  the  results  for  the  lower  density  ratios  to  be  applied  for  density  ratios  up  to  about  2.  The 
evolution  has  been  examined  for  several  Reynolds  and  Weber  numbers. 
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A.  Results 


Figure  20  shows  the  evolution  of  drops  with  We  =  2.73, 13.7, 27.4, 54.7,  and  oo.  The  density 
ratio  is  pd/p0  =  1.15,  the  Reynolds  number  is  fixed  at  331,  and  the  viscosity  ratio  is  1.  The 
computation  was  done  using  a  128  x  256  grid  and  the  drop  contours  are  shown  at  selected  times. 
The  low  Weber  number  drop  in  (a)  oscillates  due  to  the  high  surface  tension.  As  We  is  increased 
to  13.7,  the  drop  starts  to  develop  an  indentation  at  the  top  as  shown  in  (b)  but  as  it  falls,  the 
momentum  of  the  drop  decreases  and  the  surface  tension  causes  it  to  oscillate.  The  drop  shown 
in  (c)  for  We  =  27.4  also  deforms  into  an  indented  ellipsoid  initially.  The  indentation  deepens 
progressively  and  later  meets  the  bottom  of  the  drop  interface,  forming  a  forward-facing  bag  as 
observed  in  the  continuous  acceleration  case.  However,  since  there  is  no  force  to  maintain  the 
motion,  the  drop  eventually  resumes  its  initial  shape.  The  drop  with  We  =  54.7  in  (d),  on  the 
other  hand,  shows  increased  initial  deformation  which  results  in  a  bigger  rim  that  is  connected 
by  a  forward-facing  bag.  In  case  of  zero  surface  tension  in  (e),  the  drop  displays  a  roll-up  of  the 
interface.  A  similar  roll-up  has  been  observed  in  the  constant  acceleration  cases  for  drops  with  no 
surface  tension. 

In  Figure  21,  the  normalized  aspect  ratio,  centroid  velocity,  and  surface  area  are  plotted  versus 
t*  in  (a)-(c),  respectively  for  the  drops  shown  in  Figure  20.  The  aspect  ratio  plot  shows  shape  os¬ 
cillation  for  the  two  lower  values  of  We.  For  higher  W e,  the  aspect  ratio  decreases  monotonically 
to  zero  as  the  indentation  deepens  progressively.  The  velocity  plot  in  (b)  shows  a  rapid  decline 
initially.  Later,  the  velocity  of  the  drops  with  We  —  2.73  and  13.7  decreases  but  with  fluctua¬ 
tions.  The  velocities  of  the  drops  with  three  higher  Weber  numbers  decrease  monotonically  until 
it  reaches  a  minimum  and  starts  to  increase  again.  The  surface  area  plot  in  (c)  show  the  effect  of 
We  on  deformation:  As  We  increases,  the  slope  of  the  curve  increases. 

In  Figure  22,  the  evolution  of  a  drop  with  Pd/Po  =  10  is  shown  for  Re  =  242  and  pd/^o  = 
1.25.  Results  for  ten  We  =  3.74, 12.5, 18.7, 28.1, 37.4, 46.8, 56.1, 74.8, 93.5  and  oo  are  compared. 
The  computation  was  done  using  a  256  x  512  grid.  The  effect  of  increasing  We  is  generally  similar 
to  the  small  density  ratio  case.  The  drop  in  (a)  with  We  =  3.74  shows  oscillatory  deformation. 
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The  drops  with  We  =  12.5  and  18.7  in  (b)  and  (c)  develop  an  indentation  at  the  top  and  deform 
into  a  forward-facing  bag  shape.  As  the  velocity  continues  to  decrease,  the  deformation  of  the 
drops  eventually  is  reduced.  At  later  stages  of  the  deformation,  the  drops  assume  a  shape  concave 
to  the  incoming  flow  or  a  backward-facing  bag.  The  Weber  numbers  for  the  drops  are  in  the 
range  where  a  bag  breakup  mode  is  observed  experimentally  for  higher  density  ratios.  In  order 
for  the  backward-facing  bag  to  grow,  higher  initial  momentum  is  necessary.  The  drop  in  (d)  with 
We  =  28.1  shows  more  clearly  the  development  of  a  backward-facing  bag.  The  results  for  even 
higher  Weber  numbers  shown  in  (e)-(i)  show,  on  the  other  hand,  the  formation  of  a  forward-facing 
bag.  The  drop  with  zero  surface  tension  in  (j)  also  displays  a  forward-facing  bag.  In  this  case 
small  scale  irregularities  are  observed  both  on  the  edge  and  on  the  top  surface  of  the  drop. 

The  aspect  ratio,  centroid  velocity  and  surface  area  of  some  of  the  drops  shown  in  Figure  22 
are  plotted  versus  t*  in  Figure  23.  The  aspect  ratio  plot  displays  an  oscillation  of  the  drop  with 
the  lowest  We  and  monotonic  decrease  for  the  other  drops.  When  the  drops  with  We  —  12.5 
and  18.7  start  to  resume  their  initial  shape,  an  abrupt  increase  in  the  aspect  ratio  is  observed.  The 
velocity  plot  shows  a  monotonic  decline  for  all  We  shown.  Unlike  the  result  for  the  small  density 
ratio  shown  in  Figure  21(b),  no  fluctuations  occurs.  The  surface  area  plot  shows  that  the  rate  of 
deformation  increases  with  We. 

In  Figure  24  and  25,  vorticity  contours  (left)  and  streamlines  (right)  at  selected  times  are  plotted 
along  with  the  drop  contour  for  the  We  =  28.1  and  We  =  93.5  drops  shown  in  Figure  22(d)  and 
(i)  respectively.  In  both  cases,  the  vorticity  plots  show  that  most  of  the  vorticity  is  created  at 
the  outer  edge  of  the  drop,  as  expected.  The  streamline  pattern  for  the  W e  =  28.1  drop  shows 
that  the  backward-facing  bag  is  stretched  upward  in  the  downstream  direction  when  the  wake  in 
the  downstream  detaches  from  the  drop.  On  the  other  hand,  the  closed  streamlines  around  the 
We  =  93.5  drop  suggests  that  the  drop  moves  as  a  vortex  ring,  forming  a  forward-facing  bag. 

In  order  to  see  the  effect  of  the  viscosity  of  the  surrounding  fluid,  another  computation  has  been 
done  for  Re  =  387.  The  results  are  presented  in  Figure  26  for  the  same  set  of  Weber  numbers 
(We  =  3.74, 12.5, 18.7, 28.1, 37.4, 46.8, 56.1, 74.8, 93.5,  and  oo)  as  in  Figure  22.  pd/p0  is  10  and 
pd/p0  is  2.  Despite  the  change  in  the  Reynolds  number  and  the  viscous  ratio,  the  overall  evolution 
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of  the  drops  is  very  similar.  The  only  difference  is  in  the  small  structure  of  the  drop  contour.  The 
aspect  ratio,  centroid  velocity,  and  surface  area  plotted  versus  t*  for  selected  cases  in  Figure  27 
also  display  very  similar  behavior  as  in  Figure  23. 

Two  other  series  of  computations  are  shown  for  Reynolds  numbers,  Re  =  121  in  Figure  28 
and  Re  =  60.5  in  Figure  29.  The  viscosity  ratios,  Pd/Po,  are  0.625  and  0.3125,  respectively.  The , 
density  ratio  is  10  in  both  cases.  Again,  results  are  presented  for  the  same  set  of  Weber  numbers 
as  used  in  Figures  22  and  26.  Comparing  the  results  in  Figures  22, 26,  28,  and  29,  it  is  clear  that 
progressively  higher  Weber  numbers  are  necessary  in  order  to  observe  the  same  mode  of  defor¬ 
mation,  as  the  Reynolds  number  decreases.  The  translation  of  the  boundaries  between  different 
deformation  modes — oscillation,  backward-facing  bag,  and  forward-facing  bag — to  higher  We  is 
clearly  due  to  the  increased  viscous  dissipation. 

In  Figure  30,  the  effect  of  drop  viscosity  is  shown  for  drops  with  pd/po  —  10  and  Re  =  242. 
The  figures  in  each  frame  denotes  the  dimensionless  time  when  the  drop  is  plotted.  In  the  top  row, 
four  cases  are  compared  for  different  drop  viscosity  (represented  by  the  Reynolds  number  based  on 
drop  properties)  at  a  fixed  Weber  number,  We  —  28.1.  As  the  drop  viscosity  is  increased  from  left 
to  right  by  an  order  of  103,  the  drop  deformation  is  greatly  reduced.  The  least  viscous  drop  with 
Red  =  1.935  x  103  clearly  shows  the  formation  of  a  backward-facing  bag,  while  the  most  viscous 
drop  with  Red  =  1.935  remains  in  an  oblate  shape.  In  the  middle  row,  a  similar  comparison  is 
made  for  We  =  56.1.  Again,  by  increasing  the  drop  viscosity,  the  drop  deformation  is  reduced 
and  the  drop  changes  from  a  forward-facing  bag  to  an  oblate  drop.  The  result  for  We  =  93.5 
drops  shown  in  the  bottom  row  also  displays  a  similar  trend. 

Figure  31  presents  another  study  of  the  effect  of  viscosity  ratio  at  Re  =  387  and  pd/p0  =  10. 
The  overall  trend  is  similar  to  that  observed  in  Figure  30  for  Re  =  242.  As  the  drop  viscosity 
is  increased  at  a  fixed  Weber  number,  the  large  deformation  in  either  backward-facing  bag  or 
forward-facing  bag  mode  disappears. 

In  Figure  32,  a  breakup  map  is  shown  to  summarize  the  deformation  of  drops  with  pd/ p0  =  10 
and  Red  =  1935.  The  horizontal  and  vertical  axes  represent  the  Reynolds  number  and  the  Weber 
number  based  on  the  ambient  fluid,  respectively.  The  various  breakup  modes  are  denoted  by 
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different  symbols.  When  the  Weber  number  is  low,  surface  tension  prevents  large  deformation  so 
the  drop  only  oscillates.  As  the  Weber  number  increases  past  a  critical  value  (approximately  16),  a 
backward-facing  bag  starts  to  emerge  at  higher  Reynolds  numbers  first.  When  the  Weber  number  is 
higher  than  approximately  30,  the  drops  break  up  in  a  backward-facing  bag  mode  for  low  Reynolds 
numbers  and  in  a  forward-facing  bag  mode  at  high  Reynolds  numbers.  The  transitions  from  the 
backward-facing  bag  mode  to  the  forward-facing  bag  mode  occur  at  progressively  lower  Reynolds 
numbers  as  the  Weber  number  increases.  This  trend  continues  until  the  Weber  number  reaches 
approximately  100  and  all  drops  break  up  in  the  forward-facing  bag  mode  for  the  Reynolds  number 
range  examined.  Finally,  when  the  surface  tension  vanishes  (W e  — >•  oo),  the  strong  shear  due  to 
the  outside  flow  peels  off  the  drop  interface  and  shear  breakup  is  observed  when  the  Reynolds 
number  is  greater  than  100. 


B.  Conclusion 

To  study  the  characteristics  of  impulsively  accelerated  drops,  numerical  simulations  have  been 
done  for  two  density  ratios,  1.15  and  10.  These  values  of  pd/p0  are  lower  than  those  used  in  most 
of  the  experimental  investigations  of  the  drop  breakup  due  to  impulsive  disturbance  and  therefore 
of  more  relevance  to  high  pressure  sprays.  At  low  We  (<  10),  the  drops  display  oscillatory  defor¬ 
mation.  As  We  increases,  an  indentation  develops  at  the  top  of  the  drops.  Since  the  velocity  and 
the  aerodynamic  forces  causing  deformation  continue  to  decrease,  the  surface  tension  eventually 
takes  over  and  the  drops  resume  the  initial  spherical  shape. 

The  formation  of  a  backward-facing  bag  is  observed  only  for  a  finite  density  ratio  (10).  The 
absence  of  a  growing  backward-facing  bag  confirms  the  general  observation  that  the  disruptive 
aerodynamic  force  must  be  imposed  for  a  sufficiently  long  duration  of  time,  as  in  the  case  of 
a  continuously  accelerating  drop  (for  example,  see  Sadhal,  Ayyaswamy,  and  Chung28).  In  our 
simulations,  the  velocity  decreases  too  fast  for  the  backward  facing  bag  to  grow  when  the  density 
ratio  is  low. 

When  We  is  further  increased,  the  initial  deformation  is  so  large  that  a  forward-facing  bag  is 
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formed.  It  is  expected  that  in  order  to  obtain  a  shear  breakup  where  a  film  is  stripped  from  the 
main  drop,  simulations  for  even  higher  density  ratio  is  necessary. 


V.  SHEAR  BREAKUP  OF  IMMISCIBLE  FLUID  INTERFACES 

To  start  looking  at  the  shear  breakup  of  jets,  several  two-dimensional  simulations  of  immis¬ 
cible  periodic  shear  layers  were  done.  The  Reynolds  numbers  were  selected  to  be  sufficiently 
high  so  that  the  initial  instability  was  well  predicted  by  inviscid  theory,  but  viscous  effects  became 
important  at  larger  amplitude.  The  linear  stability  analysis  predicts  that  surface  tension  stabilizes 
short  waves  and  yields  a  wavelength  with  a  largest  linear  growth  rate  (most  unstable  wave).  Gen¬ 
erally,  it  was  found  that  the  inviscidly  most  unstable  mode  saturates  quickly  and  perturbations  of 
longer  wavelength  are  the  ones  that  grow  to  larger  amplitude.  Exactly  which  wavelength  is  the 
most  dangerous  one  depends  on  the  Reynolds  number.  Two  sets  of  simulations  were  conducted, 
one  at  zero  density  differences  and  the  other  at  density  ratios  of  ten.  For  zero  stratification,  sur¬ 
face  tension  prevents  Kelvin-Helmholtz  roll-up  as  seen  for  miscible  fluids  and  fingers  of  one  fluid 
penetrate  the  other  fluid.  The  slope  of  these  fingers  depends  on  the  nondimensional  wavelength 
(Weber  number).  While  viscous  effects  limit  the  growth  of  the  fingers  at  low  Weber  numbers,  high 
Weber  number  fingers  can  become  very  long.  At  even  higher  Weber  numbers  the  interface  start  to 
exhibit  a  behavior  more  similar  to  the  classical  nonlinear  Kelvin-Helmholtz  instability  and  rollup. 
The  transition,  however,  is  complex  and  intermediate  states  where  the  interface  folds  over  once 
before  being  stretched  into  a  long  finger  were  found,  for  example.  At  larger  density  ratios,  the 
evolution  is  no  longer  symmetric  and  waves  of  the  heavy  liquid  grow  into  the  lighter  one.  As  for 
zero  stratification,  waves  with  wavelength  close  to  the  most  unstable  one  are  generally  stabilized 
at  large  amplitude  by  viscous  effects  but  longer  wavelengths  lead  to  a  “wave  breaking”  where  a 
finger  of  the  heavy  fluid  is  pulled  into  the  lighter  fluid.  Even  in  the  two-dimensional  simulations, 
these  fingers  eventually  break  down  into  drops.  The  slope  of  the  “fingers”  and  the  size  of  the 
resulting  drops  depend  on  the  nondimensional  wavelength  (Weber  number). 

These  preliminary  studies  have  been  described  in  G.  Tryggvason  and  S.O.  Unverdi.  ‘The 
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Shear  Breakup  of  an  Immiscible  Fluid  Interface”  published  in  “Fluid  Dynamics  at  Interfaces”  (Ed. 
W.  Shyy)  Cambridge  University  Press,  1999.  A  more  detailed  examination  of  the  breakup  of 
immiscible  fluid  interfaces  has  been  conducted  under  an  AASERT  contract  (F49620-97-1-0525). 

VL  THREE-DIMENSIONAL  SIMULATIONS  OF  DROP  BREAKUP 

While  axisymmetric  simulations  capture  well  the  initial  evolution  of  drops  that  are  breaking  up, 
the  eventual  breakup  is  a  fully  three-dimensional  process.  During  bag  breakup,  the  rim  generally 
undergoes  Rayleigh-Taylor  instability  and  forms  drops  that  are  considerably  larger  than  drops  that 
are  formed  when  the  bag  breaks.  The  sheet  that  is  pulled  from  the  edge  of  drops  undergoing  shear 
breakup  is  also  likely  to  become  three-dimensional  before  it  breaks  into  drops.  In  some  cases, 
particularly  at  low  Eotvos  and  Weber  numbers  and  in  the  transition  zone  between  a  bag  and  a 
shear  breakup  the  process  can  be  fully  three-dimensional  right  from  the  start.  To  begin  examining 
the  three-dimensional  aspects  of  drop  breakup,  several  computations  have  been  done  of  drops 
accelerated  by  a  constant  body  force  in  both  the  bag  and  the  shear  breakup  mode.  Drops  in  bag 
breakup  mode  become  three-dimensional  quickly,  but  so  far,  the  computations  have  not  shown  the 
development  of  three-dimensional  structures  in  the  shear  breakup  mode. 

Figure  33  shows  four  frames  from  one  simulation  of  a  drop  undergoing  bag  breakup  and  the 
development  of  three-dimensional  instabilities.  The  computation  is  done  in  a  fully  periodic  domain 
on  a  643  grid.  To  resolve  the  drop  as  well  as  possible,  it  is  taken  to  be  0.4  times  the  size  of  the 
domain.  The  density  ratio  here  is  2,  and  the  Ohnesorge  numbers  are  Ohd  =  0.5  and  Oh0  =  0.3536. 
The  Eotvos  number  is  Eo  =  160  and  since  the  fluids  are  very  viscous  the  drop  breaks  up  in  a  bag- 
breakup  mode.  In  the  first  frame  the  drop  has  developed  an  indentation  at  the  back  that  grows 
until  the  drop  consists  of  a  “bowl”  with  a  thin  bottom  and  a  thick  rim.  The  bottom  of  the  “bowl 
is  pushed  back  by  the  difference  between  the  stagnation  pressure  at  the  front  of  the  drop  and  the 
low  pressure  in  the  wake.  As  the  bag  expands,  the  rim  becomes  unstable  and  starts  breaking  up 
into  few  relatively  large  drops.  The  formation  of  small  scales  eventually  results  in  structures  that 
are  not  well  resolved,  and  the  computations  are  terminated  when  this  happens.  To  compute  the 
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breakup  for  lower  Ohnesorge  numbers,  finer  resolution  is  necessary. 


VII.  SIMULATIONS  OF  HEAT  TRANSFER  DURING  BREAKUP 


The  ultimate  reason  for  atomizing  a  fuel  jet  is  to  increase  the  evaporation  rate  of  the  fuel. 
Before  evaporation,  the  fuel  must  heat  up.  As  fuel  drops  break  up,  their  surface  area  increases 
greatly  and  therefore  the  heat  transfer  to  the  drop.  To  examine  the  heat  transfer,  several  simulations 
have  been  done  that  include  the  solution  of  the  unsteady  energy  equation.  If  heat  generation  due 
to  viscous  friction  and  radiation  heat  transfer  is  neglected,  the  energy  equation  is: 

^^  +  V-(pcpT)  =  V-kVT.  (26) 


This  equation  is  solved  on  a  fixed  grid  by  an  explicit  second  order  method  in  the  same  way  as 
the  momentum  equation.  Figure  34(a)  shows  average  temperature  of  the  three  drops  versus  time 
and  Figure  34(b)  shows  the  normalized  heat  transfer  rate,  obtained  by  differentiating  the  data  in 
Figure  34(a)  and  dividing  by  temperature: 


1  dT 
h~Toa-T  dt 


(27) 


The  Ohnesorge  number  is  0.05,  Pd/ p0  =  1.15,  Pd/ Po  =  0.2,  cVd/cPo  =  1,  kd/k0  =  0.2,  and  the 
Eotvos  number  is  2.4,  24,  and  144.  For  the  lowest  Eo,  the  drop  remains  almost  spherical,  for 
Eo  =  24  the  drop  forms  a  backward  facing  bag,  and  for  the  highest  Eo,  the  drop  evolves  in  the 
shear  breakup  mode.  Since  the  highest  Eo  drop  has  reached  the  temperature  of  the  ambient  fluid 
by  time  20,  the  heat  transfer  coefficient  has  not  been  computed  after  that  time.  Initially  the  drops 
gain  heat  at  approximately  the  same  rate,  but  as  they  start  to  deform,  the  rate  depends  strongly 
on  the  breakup  mode.  The  heat  transfer  rate  is  slowest  for  the  nearly  spherical  drop  and  quickly 
increases  as  the  Eo  increases.  For  the  spherical  drop  the  heat  transfer  rate  decreases  as  the  thermal 
boundary  layer  grows  and  then  oscillates  weakly  as  the  drop  shape  oscillates.  The  heat  transfer 
coefficient  for  the  drop  that  evolves  into  a  bag-breakup  mode  first  decreases  and  then  increases 
slightly.  The  heat  transfer  for  the  higher  Eo  drops  initially  increases  rapidly  as  the  drop  deforms, 
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and  even  though  it  drops  again,  it  still  remains  about  three  times  higher  than  for  the  lower  Eo 
drops. 

The  temperature  field  and  the  drop  contour  are  plotted  for  all  three  cases  in  Figure  35.  The 
spherical  drop  is  plotted  at  t  =  28,  the  bag-breakup  drop  at  t  =  28,  and  the  shear  breakup  drop  at 
t  =  4.  The  temperature  field  around  the  spherical  drops  shows  that  recirculation  inside  the  drops 
leads  to  a  ring  of  cold  fluid  near  the  front  of  the  drop  and  a  steep  temperature  gradient  near  the 
front  of  the  drop.  As  the  drop  heats  up,  the  ambient  fluid  cools  down  and  a  cold  thermal  wake 
extends  downstream  from  the  drop.  The  temperature  field  inside  the  drop  undergoing  bag  breakup 
is  very  different.  There  is  essentially  no  recirculation  inside  the  drop  and  the  coldest  fluid  is  at  the 
centerline  of  the  drop.  Since  there  is  a  large  wake  behind  the  drop,  the  cold  thermal  wake  is  much 
larger  than  for  the  low  Eo  drop.  The  drop  undergoing  shear  breakup  has  only  a  small  thermal 
wake,  but  the  large  deformation  of  the  drop  leads  to  a  large  heat  transfer  rate.  As  the  results  in 
section  C  showed  the  rate  of  increase  of  surface  area  increases  greatly  as  Eo  is  increased.  However, 
the  flow  field  is  also  different  and  convective  heat  transfer  changes.  How  much  of  the  increased 
heat  transfer  can  be  attributed  to  increased  surface  area  and  how  much  to  increased  convection  has 
not  been  examined  in  detail  yet. 
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FIGURES 


FIG.  1.  Schematic  illustration  of  the  computational  setup. 

FIG.  2.  Resolution  test.  The  breakup  of  a  drop  computed  using  a  256  x  512  grid  (left)  and  a  512  x  1024 
grid  (right).  pdlpo  =  1.15,  Eo  -  144,  Oh0  =  0.05,  Ohd  =  0.0466.  The  drop  shape  is  plotted  every 
At*p  =  3.873. 

FIG.  3.  Resolution  test.  Aspect  ratio  and  centroid  velocity  plotted  versus  t* .  Results  us¬ 
ing  three  different  grids,  128  x  256,  256  x  512,  and  512  x  1024,  are  shown.  pd/p0  =  1.15, 
Eo  =  144,  Oh0  =  0.05,  Ohd  =  0.0466. 

FIG.  4.  Test  of  the  Boussinesq  approximation.  Aspect  ratio  and  centroid  velocity  plotted  versus  t. 
Results  are  shown  for  four  different  Atwood  numbers:  0.07,  0.11,  0.2,  and  0.33.  The  corresponding  density 
ratios  are  1.15,  1.25, 1.5,  and  2.0.  Eo  =  72,  Oh0  =  0.241,  and  pd/p0  =  1. 

FIG.  5.  Effect  of  Eo  on  the  deformation  of  drops  with  pd/ p0  =  10.  Oh0  =  Ohd  =  0.05.  The  simula¬ 
tions  are  done  using  a  256  x  768  grid  for  a  moving  computational  domain  of  dimensions  5  x  15  the  initial 
drop  diameter.  The  boundaries  of  the  column  do  not  indicate  the  actual  boundaries  of  the  computational 
domain.  The  gap  between  two  successive  drops  in  each  column  represents  the  distance  the  drop  travels  at 
a  fixed  time  interval  and  the  last  interface  is  plotted  at  t*  =  (a)  11.19;  (b)  15.82;  (c)  14.85;  (d)  13.83;  (e) 
11.19;  (f)  7.15;  (g)  7.83;  (h)  6.78;  (i)  5.54. 

FIG.  6.  Effect  of  Eo  on  the  deformation  of  drop  with  pd/ p0  =  1.15.  OhQ  =  0.05,  Ohd  —  0.0466.  The 
simulations  are  done  using  a  256  x  1280  grid  in  (b)  and  (c)  and  a  256  x  768  grid  in  (a),  (d),  and  (e).  The 
fixed  computational  domain  has  a  dimension  of  5  x  15  the  initial  drop  diameter  in  (b)  and  (c)  and  5  x  25 
the  initial  drop  diameter  in  (a),  (d),  and  (e).  The  dashed  line  in  (a),  (d),  and  (e)  represents  the  actual  bottom 
boundary  of  the  computational  domain.  The  gap  between  two  successive  drops  in  each  column  represents 
the  distance  the  drop  travels  at  a  fixed  time  interval  and  the  last  interface  is  plotted  at  t*  =  (a)  44.72;  (b) 
79.06;  (c)  89.44;  (d)  37.94;  (e)  38.73. 
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FIG.  7.  Vorticity  contours  (left)  and  streamlines  (right)  for  the  drop  in  Figure  5(c).  pd/ Po  —  10, 
Eo  =  28.8,  Oh0  =  Ohd  =  0.05.  The  results  are  shown  for  four  selected  times. 

FIG.  8.  Vorticity  contours  (left)  and  streamlines  (right)  for  the  drop  in  Figure  5(f).  pd/ po  =  10, 
Eo  =  60,  Oh0  =  Ohd  =  0.05.  The  results  are  shown  for  five  selected  times. 

FIG.  9.  Vorticity  contours  (left)  and  streamlines  (right)  for  the  drop  in  Figure  6(e).  Pd/Po  =  1*15, 
Eo  -  144,  Oh0  =  0.05,  Ohd  =  0.0466.  The  results  are  shown  for  four  selected  times. 

FIG.  10.  Centroid  velocity  versus  t*  for  the  drops  shown  in  Figure  5.  The  results  are  presented  for 
Eo  =  12, 24, 28.8, 36, 48, 60, 72, 96,  and  144.  pd/p0  =  10,  Oh0  =  Ohd  =  0.05. 

FIG.  11.  Centroid  velocity  versus  t*  for  the  drops  shown  in  Figure  6.  The  results  are  presented  for 
Eo  =  12, 24,48, 96,  and  144.  pd/ p0  =  1.15,  Oha  =  0.05,  Ohd  =  0.0466. 

FIG.  12.  Normalized  surface  area  versus  t*  for  the  drops  shown  in  Figure  5.  The  results  are  presented 
for  Eo  =  12, 24, 28.8, 36, 48, 60, 72, 96,  and  144.  pd/ po  =  10,  Oh0  =  Ohd  =  0.05. 

FIG.  13.  Normalized  surface  area  versus  t*  for  the  drops  shown  in  Figure  6.  The  results  are  presented 
for  Eo  =  12,24,48, 96,  and  144.  pd/p0  =  1.15,  0ho  =  0.05,  Ohd  =  0.0466. 

FIG.  14.  Effect  of  Oh  for  drops  with  pd/po  =  10-  The  drop  evolution  is  shown  for  three 
Oh0  =  Ohd  =  0.05,  0.125,  and  0.25.  In  the  upper  row  (a)-(c),  Eo  is  fixed  at  28.8  and  the  time  inter¬ 
val  between  successive  interfaces,  A£*.  is  2.121.  In  the  lower  row  (d)-(f),  Eo  is  144  and  At*  is  0.791. 

FIG.  15.  Effect  of  Oh  on  the  initial  deformation  of  drops  with  Pd/Po  =  1-15.  In  all  cases,  Pd/Po  =  1 
and  the  time  intervals  between  successive  interfaces.  At*  =  1.581.  In  the  upper  row  (a)-(d),  Eo  =  24  and 
in  the  lower  row  (e)-(h),  Eo  =  oo  (zero  surface  tension).  The  viscous  simulations  (a)-(c)  and  (eMg)  were 
done  using  a  128  x  384  grid,  (d)  and  (h)  were  done  using  an  inviscid  vortex  method. 
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FIG.  16.  Effect  of  the  viscosity  ratio.  In  each  row,  the  drops  at  a  late  time  are  shown  for  different  values 
of  the  drop  viscosity  while  other  parameters  are  fixed.  In  (a)  and  (b),  Ohd  =  0.05, 0.125, 0.25,  and  1.25 
(from  left  to  right).  In  (c)  and  (d),  Ohd  =  0.0093, 0.0466, 0.2331,  and  1.1656  (from  left  to  right). 

FIG.  17.  Centroid  velocity  versus  t*  for  the  drops  shown  in  Figure  16. 

FIG.  18.  Deformation  and  breakup  regime  maps  for  pd/p0  =  10-  Three  maps  are  shown  for 
OhQ  =  0.05,  0.125  and  0.25.  In  each  map,  the  horizontal  and  vertical  axes  are  Ohd  and  Eo,  respectively. 

FIG.  19.  Deformation  and  breakup  regime  maps  for  pd/po  =  1-15.  Three  maps  are  shown  for 
Oh0  =  0.05, 0.25  and  1.25.  In  each  map,  the  horizontal  and  vertical  axes  are  Ohd  and  Eo,  respectively. 

FIG.  20.  Evolution  of  impulsively  started  drops  with  Re  =  331,  pd/ p0  =  1-15,  and  pd/p0  =  1- 
Results  for  five  We  are  shown  as  denoted  by  the  numbers  below  the  figures.  The  numbers  inside  the  frames 
denote  t*  when  the  interfaces  are  plotted.  The  computations  were  done  on  a  128  x  256  grid. 

FIG.  21.  Non-dimensionalized  aspect  ratio,  centroid  velocity,  and  surface  area  plotted  versus  t*  for  the 
drops  shown  in  Figure  20.  Results  for  five  We  =  2.73,  13.7,  27.4,  54.7,  and  oo  are  compared.  Re  =  331, 
Pd/po  =  1.15  and  pd/p0  =  1. 

FIG.  22.  Evolution  of  impulsively  started  drops  with  Re  =  242,  pd/ p0  =  10,  and  pd/p0  =  1-25. 
Results  for  ten  We  are  shown  as  denoted  by  the  numbers  below  the  figures.  The  numbers  next  to  the  drops 
denote  t*  when  the  interfaces  are  plotted.  The  centroid  of  the  drops  in  a  column  are  separated  by  a  fixed 
distance.  The  gap  between  two  successive  drops  in  a  column  does  not  represents  the  distance  the  drop 
travels  during  the  time  interval.  The  computations  were  done  on  a  256  x  512  grid. 

FIG.  23.  Non-dimensionalized  aspect  ratio,  centroid  velocity,  and  surface  area  plotted  versus  t*  for 
selected  cases  of  the  drops  shown  in  Figure  22.  Results  for  five  We  =  3.74,  12.5,  18.7,  37.2,  and  93.5  are 
compared.  Re  =  242,  pd/ pQ  =  10  and  pd/Po  =  1-25. 


FIG.  24.  Vorticity  contours  (left)  and  streamlines  (right)  for  the  drop  shown  in  Figure  22(d). 
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FIG.  25.  Vorticity  contours  (left)  and  streamlines  (right)  for  the  drop  shown  in  Figure  22(i). 

FIG.  26.  Evolution  of  impulsively  started  drops  with  Re  —  387,  Pd/Po  =  10.  and  Pd/Po  =  2.  Results 
for  ten  We  are  shown  as  denoted  by  the  numbers  below  the  figures.  The  numbers  next  to  the  drops  denote 
t*  when  the  interfaces  are  plotted.  The  centroid  of  the  drops  in  a  column  are  separated  by  a  fixed  distance. 
The  gap  between  two  successive  drops  in  a  column  does  not  represents  the  distance  the  drop  travels  during 
the  time  interval.  The  computations  were  done  on  a  256  x  512  grid. 

FIG.  27.  Non-dimensionalized  aspect  ratio,  centroid  velocity,  and  surface  area  plotted  versus  t*  for 
selected  cases  of  the  drops  shown  in  Figure  26.  Results  for  five  We  =  3.74,  12.5,  18.7,  37.2,  and  93.5  are 
compared.  Re  —  387,  Pd/po  —  10  and  pd/po  —  2. 

FIG.  28.  Evolution  of  impulsively  started  drops  with  Re  =  121,  Pd/Po  =  10,  and  pd/p<>  —  0.625. 
Results  for  ten  We  are  shown  as  denoted  by  the  numbers  below  the  figures.  The  numbers  next  to  the  drops 
denote  t*  when  the  interfaces  are  plotted.  The  centroid  of  the  drops  in  a  column  are  separated  by  a  fixed 
distance.  The  gap  between  two  successive  drops  in  a  column  does  not  represents  the  distance  the  drop 
travels  during  the  time  interval.  The  computations  were  done  on  a  256  x  512  grid. 

FIG.  29.  Evolution  of  impulsively  started  drops  with  Re  =  60.5,  Pd/Po  =  10,  and  pd/po  =  0.3125. 
Results  for  ten  We  are  shown  as  denoted  by  the  numbers  below  the  figures.  The  numbers  next  to  the  drops 
denote  t*  when  the  interfaces  are  plotted.  The  centroid  of  the  drops  in  a  column  are  separated  by  a  fixed 
distance.  The  gap  between  two  successive  drops  in  a  column  does  not  represents  the  distance  the  drop 
travels  during  the  time  interval.  The  computations  were  done  on  a  256  x  512  grid. 

FIG.  30.  Drop  viscosity  effect  on  the  deformation  of  drops  with  pd/ Po  —  10  and  Re  =  242.  The 
number  in  each  frame  denotes  the  dimensionless  time  when  the  drop  is  plotted. 

FIG.  31.  Drop  viscosity  effect  on  the  deformation  of  drops  with  pd/Po  —  10  and  Re  =  387.  The 
number  in  each  frame  denotes  the  dimensionless  time  when  the  drop  is  plotted. 

FIG.  32.  Breakup  mode  map  for  impulsively  started  drops  with  Pd/Po  =  10  and  Red  =  1935. 
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FIG.  33.  Three-dimensional  simulation  of  drop  breakup,  pd/ Po  =  2,  Eo  =  160,  Ohd  —  0.5,  and 
Oh0  =  0.3536. 

FIG.  34.  Simulations  of  heat  transfer  during  drop  breakup  for  Eo  =  2.4, 24,  and  144.  Average  tempera¬ 
ture  of  the  drops  (a)  and  normalized  heat  transfer  rate  (b)  are  plotted  versus  time.  Oh  =  0.05,  Pd/ Po  =  1-15, 
Pdf  Po  ~  0-2,  cvJcVo  =  1,  and  kd/k0  0.2 

FIG.  35.  Temperature  field  and  drop  contour  are  plotted  for  the  drops  shown  in  Figure  34  at  selected 
times. 
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